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ABSTRACT 

We investigate numerically the role of thermal instability (TI) as a generator of 
density structures in the interstellar medium (ISM), both by itself and in the context of 
a globally turbulent medium. We consider three sets of numerical simulations: a) flows in 
the presence of the instability only; b) flows in the presence of the instability and various 
types of turbulent energy injection (forcing), and c) models of the ISM including the 
magnetic field, the Coriolis force, self-gravity and stellar energy injection. Simulations in 
the first group show that the condenstion process which forms a dense phase ( "clouds" ) 
is highly dynamical, and that the boundaries of the clouds are accretion shocks, rather 
than static density discontinuities. The density histograms (PDFs) of these runs exhibit 
either bimodal shapes or a single peak at low densities plus a slope change at high 
densities. Final static situations may be established, but the equilibrium is very fragile: 
small density fluctuations in the warm phase require large variations in that of the cold 
phase, probably inducing shocks into the clouds. Combined with the likely disruption 
of the clouds by Kelvin- Helmholtz instability (Murray et al. 1993), this result suggests 
that such configurations are highly unlikely. 

Simulations in the second group show that large-scale turbulent forcing is incapable 
of erasing the signature of the TI in the density PDFs, but small-scale, stellar-like 
forcing causes the PDFs to transit from bimodal to a single-slope power law, erasing the 
signature of the instability. However, these simulations do not reach stationary regimes, 
the TI driving an ever-increasing star formation rate. Simulations in the third group 
show no significant difference between the PDFs of stable and unstable cases, and reach 
stationary regimes, suggesting that the combination of the stellar forcing and the extra 
effective pressure provided by the magnetic field and the Coriolis force overwhelm the 
TI as a density-structure generator in the ISM, the TI becoming a second-order effect. 
We emphasize that a multi-modal temperature PDF is not necessarily an indication of 
a multi-phase medium, which must contain clearly distinct thermal equilibrium phases, 
and that this "multi-phase" terminology is often inappropriately used. 
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1. Introduction 

The idea tliat tliermal instability (TI) plays a fun- 
damental role in controlling important aspects of the 
interstellar gas and star formation in galaxies has 
been tremendously influential. As a model for the 
formation of cool dense clouds in pressure equilibrium 
with a warm rarefied intercloud medium, TI has been 
invoked to explain the existence of diffuse interstellar 
clouds (Field, Goldsmith, & Habing 1969, FGH), as 
a means of regulating the mass flow between differ- 
ent components of the ISM and the star formation 
rate (Chieze 1987; Parravano 1989), as a major fac- 
tor allowing star formation to occur at an apprecia- 
ble rate in galaxies (Norman & Spaans 1997; Spaans 
& Norman 1997), the fragmentation of protoglobu- 
lar clusters (Murray & Lin 1996), and cooling flows 
around luminous cluster elliptical galaxies and more 
general situations (see Nulsen & Fabian 1997; Allen 
k Fabian 1998; Mathews & Brighenti 2000 for recent 
discussions) . 

The linear stability analysis for TI was first worked 
out in detail by Field (1965), who clearly delineated 
the isobaric, isochoric, and iscntropic modes, and ex- 
amined effects such as magnetic fields and conduction. 
Subsequent work generalized the anlaysis to include 
chemical reactions (Yoneyama 1973), and refined the 
calculation of relevant heating and cooling rates (see 
Wolfire et al. 1995 for a recent discussion). In its 
simplest form, the isobaric mode of TI occurs when a 
compression leads to such enhanced cooling that pres- 
sure in the compressed region decreases. This leads 
to a picture of dense clouds in pressure equilibrium 
with their surroundings, the basis for the "two-phase" 
ISM model of FGH. This conceptual framework was 
extended by McKee & Ostrikcr (1977) to include a 
third, hot, phase representing the interiors of expand- 
ing supernova remnants, but still assumes that clouds 
form by TI. 

The idea of cloud formation by TI suffers from a 
serious defect: the linear stability analysis assumes a 
static uniform initial medium. It is now observation- 
ally evident that the ISM in galaxies is in a highly 
agitated, "turbulent" state, and so it is questionable 
whether the static initial conditoins are ever realized. 
The basic question is whether dynamical motions of 
the gas (shocks, cloud collisions, or the turbulence it- 
self) will disrupt incipient clouds faster than they can 
condense. 

Even within the context of linear stability analy- 



ses, there have been several suggestions that gravita- 
tionally induced convective motions could suppress TI 
(e.g. Balbus & Soker 1989). Numerical studies of the 
nonlinear evolution of TI in the solar chromosphere 
transition region and corona (Dahlberg et al. 1987; 
Karpen et al. 1988) and in galaxy cluster cooling 
flows (Loewenstein 1989; Malagoli et al. 1990; Hat- 
tori & Habe 1990; Yoshida et al. 1991) have generally 
confirmed this suspicion. The motion of the incipient 
condensations leads to vortices, or to complete dis- 
ruption by Kelvin-Helmholtz and Rayleigh- Taylor in- 
stabilities, although a magnetic field can diminish the 
effect in these convective cases. Moreover, Murray et 
al. (1993) have showed in some detail how the motion 
of a condensation initially in a two-phase medium in 
pressure equilibrium can easily lead to disruption of 
the condensation by dynamical instabilities, similar 
to the results quoted above. Strongly self-gravitating 
condensations can avoid this fate, but there is no rea- 
son to assume that such condensations were formed 
by thermal instability; supersonic turbulence interac- 
tions are sufficient to produce such structures (e.g., 
Vazquez-Scmadcni et al. 1996 ). 

What about thermal instability in a supersonic tur- 
bulent medium, like the gas in galaxies? (By "turbu- 
lent" we mean a disordered but not completely ran- 
dom velocity [and density] field covering a large range 
of scales.) In the ISM of the Milky Way there is cer- 
tainly strong evidence for supersonic motions at all 
scales above at least 0.1 pc, in every sort of envi- 
ronment. Supersonic spectral linewidths are observed 
even in regions with no detectable internal star for- 
mation (e.g. the Maddalena molecular cloud complex; 
see Williams & Blitz 1998) and in diffuse mostly-H I 
clouds in which self-gravity is unimportant (known 
since the 1950s; see Heithausen 1996 for a recent 
study of a subclass of such clouds), as well as in the 
more intensively-studied star-forming molecular cloud 
structures. 

A simple argument might suggest that thermal 
instability may be inoperative, or at least less efh- 
cient, as a cloud formation process in a supersonically- 
turbulent ISM. The isobaric mode of TI in a region 
of size L condenses on a characteristic time scale 
L/c, where c is the sound speed, while the turbu- 
lence shears (or expands, or compresses) the forming 
condensations on a crossing time scale L/v, where v 
is the characteristic turbulent speed at scale L. If the 
turbulence is supersonic, the incipient condensation 
should be disrupted faster than it can condense. How- 
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ever this argument only applies in an average sense; a 
region of incipient condensation might avoid disrup- 
tive turbulent interactions for an unusually long time. 
More importantly, as we show in sec. 3.2, when one 
considers that the instability growth time decreases 
with scale size while the crossing time increases with 
scale, the above argument can only apply at suffi- 
ciently large scales, and below some size scale the in- 
stability could still operate. Moreover, Hennebelle & 
Perault (1999) have recently shown that TI can be 
triggered by a dynamical compression in an originally 
stable flow. 

A quantitative result suggesting that the idea of 
pressure-confined clouds, central to models of the 
ISM which rely on thermal instability, as well as 
other quasi-static conceptions, is not viable in a su- 
personically turbulent medium was given, using nu- 
merical simulations, by Ballesteros-Paredes, Vazquez- 
Semadeni & Scalo (1999, herafter BVS99). They 
showed, partly through the evaluation of volume and 
surface terms in the virial theorem, that the impor- 
tance of the kinetic energy surface terms implies that 
clouds cannot be considered as quasi-permanent enti- 
ties with real "boundaries," but are instead continual- 
lly changing, forming and dissolving. Not only is ther- 
mal pressure incapable of confining turbulent den- 
sity fluctuations, but the idea of external turbulent- 
pressure confinement seems self-defeating, since the 
external turbulent stresses mostly serve to distort and 
disrupt clouds. 

There are additional problems with the two- or 
three-phase ISM models. One involves time scales. 
The time scale on which a given parcel of the ISM 
is subjected to impulsive perturbations, whether by 
heating or by shock interactions, may be smaller than 
the time scale for TI, and than the time for read- 
justment to pressure equilibrium. That clouds must 
be subjected to frequent stochastic heating events 
was the basis for Kahn's (1955) early cloud-cloud 
collision model for the temperature distribution of 
clouds, and the "time-dependent ISM" models pro- 
posed by Bottcher et al. (1970) and examined in de- 
tail by Gerola, Kafatos, & McCray (1974). The latter 
work only included the influence of stochastic heat- 
ing and ionization sources (no hydrodynamics), but 
the resulting statistical description of the diffuse ISM 
(Gerola et al. 1974) was in better agreement with 
available observations than the static two-phase mod- 
els. A recent summary of (mostly observational) argu- 
ments against the specific three-phase McKee & Os- 



triker model is given by Elmegreen (1997). Moreover, 
the conclusion that pressure equilibrium is unlikely 
because of the frequency of either cloud collisions or 
shocks from supernovae or cluster superbubbles has 
been found independently several times in the liter- 
ature (e.g. Stone 1970; Heathcote & Brand 1984; 
Bowyer et al. 1995; Berghofer et al. 1998; Kornreich 
& Scalo 2000). 

Given all the above results, it is important to 
examine whether the traditional picture of a multi- 
phase ISM structured by TI should be retained as 
a useful model. By "multi-phase" we are referring 
to a medium containing various different thermody- 
namic equilibrium regimes, some of them possibly sta- 
ble and others unstable, with the requirement that a 
multi-phase medium involves fluid parcels undergoing 
a "phase transition" when transiting between these 
equilibria. The structuring of the density field by the 
multiphase nature of the medium should be reflected 
in multimodal density and temperature pdfs, with gas 
accumulating in the stable "phases" . 

However, note that the temperature field is ex- 
pected to show a bimodal PDF even in the absence 
of TI, simply because of the functional shape of the 
interstellar cooling curve, which has the form of two 
plateaus separated by a sharp discontinuity at ~ 10^ 
K (see Dalgarno & McCray 1972, Fig. 11). For ex- 
ample, the simulations of Korpi et al. (1999) do not 
contain real phases, since no background heating ca- 
pable of balancing the cooling is included; yet, the 
temperature still exhibits a bimodal distribution in 
those simulations. But for a Tl-structured multi- 
phase ISM, the gas must be concentrated into two 
(or more) density components: dense clouds and a 
rarified intercloud medium. Two-dimensional simula- 
tions of the ISM including a typical cooling function, 
but containing no TI (Scalo et al. 1998) showed no 
indication of a bimodal density PDF. For this reason 
in this paper we consider a bimodal density distribu- 
tion as the signature of a Tl-structured ISM. Actu- 
ally, our experiments have shown that, under certain 
conditions, the PDFs of unstable flows may be uni- 
modal, but with an added slope change where the 
second peak should be located. Thus, we consider ei- 
ther of these PDF shapes as a signature of density 
structuring by TI. 

It is also of great interest to investigate whether 
TI can significantly enhance the formation of clouds 
by turbulent interactions. As shown in Vazquez- 
Semadeni, Ballesteros-Paredes & Rodriguez (1997), 
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BVS99 and references therein, "clouds" naturally 
form in a supersonic turbulent medium, without any 
need for instabilities, as turbulent velocity streams 
interact to form compressed layers, as suggested by 
Elmegreen (1993). 

The present work attempts to investigate these 
ideas using numerical simulations in two dimensions. 
First we simulate the evolution of an initially static 
gas whose cooling and heating functions should guar- 
antee TI, and show that, as expected, the density dis- 
tribution shows the signature of the instability. Sec- 
ond, we perform simulations of a randomly forced tur- 
bulent flow in the presence of the instability, showing 
that large-scale forcing does not seem to be able to 
prevent the development of the instability, but small- 
scale forcing does. 

Finally, a third series of simulations considers a 
supersonic turbulent gas, including self-gravity, the 
magnetic field, the Coriolis force, and energy input 
due to star formation, in stable, marginal and unsta- 
ble cases. We demonstrate that in the latter case the 
density distribution does not show the signature of 
the instability. The temperature and density are still 
anticorrelated, as expected from the cooling function, 
but a two-phase density field is not realized. This is 
ultimately due to the fact that, when the cooling time 
scale is much smaller than the hydrodynamic time 
scale, the density is "slaved" to the dynamics, not the 
temperature (cf., Vazquez-Semadeni, Passot & Pou- 
quet 1995); another way of expressing this is that the 
density field is controlled by ram pressure/advection, 
not the thermal pressure. 

2. The Model 

We use the numerical ISM model presented by Pas- 
sot, Vazquez-Semadeni & Pouquet (1995), which uses 
a single-fluid approach to describe the ISM on a 1- 
kpc^ plane on the Galactic disk, centered at the so- 
lar Galacto-centric distance. In this model, various 
physical ingredients can be included at will, such as 
self-gravity, the magnetic field, disk rotation, as well 
as model terms for the radiative cooling (A) , a diffuse 
background radiation (Fd), and energy input due to 
star formation (Fg). In this paper, we solve the equa- 
tions in logarithmic form for the density and temper- 
ature, as done in Gazol-Patiho & Passot (1999). This 
is because the shocks induced by the development 
of the instability are quite strong, and the logarith- 
mic approach avoids the appearance of negative val- 



ues of the density and temperature due to the Gibbs 
phenomenon in the vicinity of discontinuities, giving 
somewhat more robustness to the code. The equa- 
tions solved read 



dlnp 
dt 
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where p is the fluid density, u the velocity, P the 
thermal pressure, the gravitational potential, 11 the 
angular velocity due to Galactic rotation, B the mag- 
netic field and e the specific internal energy, related to 
the temperature T by e = CvT, Cv being the specific 
heat at constant pressure. In eq. (||), J = Lo/L^ is the 
Jeans number, measuring the length of the integration 
box in units of the Jeans length, and M = Mo/cq is 
the Mach number of the velocity unit with respect 
to the temperature unit, given by its corresponding 
speed of sound Cq. In eq. (||), 7 is the ratio of specific 
heats at constant pressure and volume, and kt is the 
thermal diffusivity. 

Momentum and magnetic field dissipation are in- 
cluded by a combination of a quadrilaplacian "hyper- 
viscosity" operator which confines the viscous effects 
to the smallest scales, and a second order operator, 
which filters out possible oscillations in the vicinity 
of strong shocks. The use of the hyperviscosity op- 
erator allows the use of much smaller (typically by 
a factor of 10) values of the second-order diffusivi- 
ties than otherwise. An artificial mass diffusion term 
in the continuity equation has been added in order to 
smooth out the density gradients, allowing the code to 
handle stronger shocks. Its effect is only to spread out 
and somewhat compress density peaks. The effect of 
this term has been quantified by Vazquez-Semadeni, 
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Ballesteros-Paredes & Rodriguez (1997) and BVS99. 

The physical units are Lq = 1 kpc, po = Imncm"^ , 
Uo = 11.7km s"\ to = 1.3 X 10^ yr, To = 10"* K, and 
Bo — 5/iG. Throughout the paper, we omit toh, the 
proton mass, in the specification of densities. 

The above equations are solved in two dimensions 
at a resolution of 128^ grid points (implying a spatial 
resolution of 7.8 pc) by means of a pseudo-spectral 
method, which imposes periodic boundary conditions. 
The temporal scheme is a third order Runge-Kutta 



for the non-linear terms, combined with a Crank- 



Nicholson for the linear terms. The initial conditions 
are set up in Fourier space, uncorrelated in all vari- 
ables, and characterized by a power spectrum of the 



form k'^ exp 



-2{k/ko) 



with random phases. In the 
case of the magnetic field, the initial condition in- 
volves a uniform component of strength Bi = 1.5/iG 
in the x-direction mimicking the azimuthal compo- 
nent of the field on the Galactic plane, upon which a 
fluctuating field of rms strength 5/iG is added. 

An important note is that, in spite of the usage 
of logarithmic variables and hyperviscosity, our sim- 
ulations cannot handle very strong shocks because of 
the spectral scheme used (the Gibbs phenomenon pro- 
duces artificial oscillations near discontinuities), and 
the refore the simulations must be stopped when the 



gra hients of the physical variables become too steep . 



However, they are able to reach advanced enough 
stages of the development of the instability that its 
fully dynamical character is seen, and the structures 
arising are well defined. 

We describe below only those model terms which 
are directly related with simulations presented in the 
following sections. Further details concerning the 
choice of the other terms and parameters can be found 
in Passot et al. (1995) and Vazquez-Semadeni et al. 
(1996). Table 1 presents the relevant parameters for 
the runs in this paper. 

2.1. Radiative cooling and diffuse heating 

The heating (F) and cooling (pA) rates are per unit 
mass, and thus differ somewhat from the standard 
definitions, which are per unit volume. In eq. (|^), 
the term F^ models the heating of the gas due to 
photoelectrons ejected from dust grains due to the 
background UV radiation field, or heating by low- 
energy cosmic rays, and is taken as a constant Yd — 
To everywhere and throughout the duration of the 
simulations. 



The simulations use a cooling function with piece- 
wise power-law dependence on the temperature of the 
form 



A = Q 



for Ti <T< n 



(7) 



The various runs have different values of Pi^i+i, given 
in Table 1. One set of values of these exponents, 
namely (3i2 = 1, /923 = 0.25, P34 = 4.18 and (345 = 
—0.53, was chosen as a reasonable piecewise power- 



law fit to the radiative cooling curve of Dalgarno & 
McCray (1972) for a fractional ionization of lO^** at 
T < IQ-* K (see also [Spitzcr 197^ ). For this fit, we 
have chosen the transition temperatures Ti as Ti = 10 
K, T2 = 398 K, T3 = 10* K, T4 = 10^ K and Tg = 10^ 
K. As discussed below, this particular combination of 
values of Pi^i+i and Ti give a thermally unstable range 
(isobaric mode) between T2 and T3. Other runs have 
variations of this cooling function in order to obtain 
thermal stability or marginal instability in the range 
398 — 10"* K without a significant increase of the cool- 
ing rate above 10"* K. Due to continuity reasons, the 
values of coefficients Cis+i must also be changed, and 
are also listed in Table 1. Finally, all runs share the 
same values of /345 and of Ti , T3 , T4 and , and are 
therefore not listed in Table 1. 

As discussed in previous papers 



e.g., Vazquez-I 



Semadeni et al. 1995| ; Vazquez-Semadeni et al. 1996| ) 



the thermal time scales are typically much shorter 
than the dynamical ones, implying that the turbulent 
motions are quasi-static compared to the background 
heating and cooling rates, and allowing for the estab- 
lishment of thermal equilibrium. With the adopted 
power-law behavior for these processes, the equilib- 
rium temperature and pressure are 



Pr, 



To 



97°, ..+1 



7 \C^, 



(8) 



(9) 



where the (piecewise) effective polytropic index is 
given by 7°^ i+i = 1 — l//3i,i+i. The isobaric mode 
of the TI, characterized by a decrease in the pressure 
as the density increases, develops when 7°i_i_|_i < 0, 
corresponding in this case to < 1. The lower 

pressure causes an effective "suctioning" action by the 
high-density regions. The values of the resulting ^"^^2 
and 7°23 are also given in Table 1. The value of 7''34 is 
always taken as 7^34 — 0.76, and therfore not listed in 
Table 1. Note that the flow behaves as P oc p"*" * in 
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the density interval Pi < p < pt+i, where the "tran- 
sition" densities p, are defined as the values whose 
corresponding equilibrium temperatures (eq. |^) co- 
incide with Note also that some of the runs, in 
addition to having a thermally unstable range with 
7*^23 < Oj have a marginally stable regime at low tem- 
peratures, with 7®]^2 = 0. This is the case in par- 



ticular for our fit of the Dalgarno fc McCray (1972 



cooling curves, and its consequences are discussed in 



§3.1.1 



The dispersion relations for the combined ther- 
mal-|-gravitational instability in the presence of a 
uniform field Bi have been derived for the two- 
dimensional case by Passot et al. (1995) (for the gen- 
eral 3D case, see Elmegreen 1994), and read: 



2 t2 7 2/7 
UJ^ = J — k [ ■ 



cq 



Bl) - 

2 2 



1 



(10) 

(11) 



where loi and Wt respectively denote the growth rates 
for perturbations in the longitudinal and transverse 
directions with respect to i?i , k denotes the wavenum- 
ber, J is the Jeans number defined above, and k = 
2U,[1 + [\/1){y/n){dy/d^l) is the epicyclic frequency. 
Unstable wavenumbers are those for which oj^ > 0. 
For the values of the parameters used in the unstable 
runs, typical unstable ranges are fc > 1 , implying that 
essentially all scales included in the simulations are 
unstable. Note that the thermal instability (7*^ < 0) 
has the effect of reversing the instability criterion with 
respect to the purely gravitational case, in the sense 
that now it is scales smaller than some threshold value 
that are unstable. 

Table 1 also contains information about the ther- 
mal equilibrium regime for the heating and cooling 
functions used here, in particular p2 and p3, and the 
equilibrium temperature T[{p)) at the mean density, 
{p) = 1 cm~'^. Note that the value of Tq used here 
is in general not the same as that used by Vazquez- 
Semadeni et al. (1995). We have chosen Tq in order 
to ensure that T{{p)) lies within the unstable range. 



^Note that runs with different cooUng functions must also have 
different transition densities P2 and P3. It is not possible to 
vary the cooling and heating functions in such a way as to 
keep both the transition densities and temperatures constant 
simultaneously, without varying the effective rates excessively. 



2.2. Forcing 

We consider two kinds of forcing (i.e., energy injec- 
tion) for the turbulent runs. The first is a large-scale 
forcing, accomplished through an additional accelera- 
tion term in the momentum equation (eq. |^), taken 
as a random vector field with a spectrum proportional 
to fc* below the forcing wavenumber fcfor = 4, and 
fc"** above fefor- We consider two subcases of this 
forcing, one with 100% solenoidal (rotational) com- 
ponents and the other with 50% solenoidal and 50% 
compressible components. 

The second type of forcing is applied at small 
scales, and consists of local, discrete heating sources 
turned on at a grid point x whenever p{x) > psF and 
V • u{x) < 0, where psF is a free parameter, but taken 
equal to 6(p) always. The sources stay on for a time 
interval At = 6 x lO^yr. This forcing mimics the ef- 
fect of ionization heating fro m massive stars, and acts 
on the scale of ~ 5 pixels ( Vazquez- Scmadcni et al 
1995| ; [Passot et al. 1995| ). 



2.3. The simulations 

In order to investigate the interplay between the 
TI and the turbulent dynamics, we have performed 
three sets of three simulations each. First, we have 
considered the development of the instability alone 
(runs labeled "pure" or, simply, "P" in Table 1), by 
simulating a purely hydrodynamic fiow without the 
magnetic field, rotaion or self-gravity. Only the ap- 
propriate cooling and background heating functions 
necessary for the existence of a thermally unstable 
range are included. The three runs differ in details of 
the cooling and background heating functions, chosen 
to show the variety of structures that develop and the 
effects on them of the closeness of the transition den- 
sities to the mean density. The label "P" is thus fol- 
lowed by a mnemonic indicating the value of Tq used 
in each run. 

The second set of simulations (referred to as the 
"hydrodynamic" runs, or "H" in Table 1) are similar 
to the P runs, except that they include a variety of 
forcing schemes, as described in §2.2. These are de- 



signed to test the effect of the scale and the amount 
of compressible content of the energy injection on the 
supression of the instability by the turbulence. As 
indicated in §^]^, the small-scale forcing is stellar- 
like, i.e., in the form of localized sources of heating 
which create expanding bubbles of hot gas, while the 
large-scale forcing is random, with varying ratios of 
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compressible-to-solenoidal content . 

Finally, we have performed three ISM-like simula- 
tions (labeled "ISM" in Table 1), including the mag- 
netic field, rotation and self-gravity, as well as the 
stellar-like forcing. The runs differ from each other in 
whether the cooling/heating combinations produce a 
stable, marginal or unstable regime, and are designed 
to test the effects of the presence of the TI in an ISM- 
like medium. 

For the "pure" and "hydro" runs we have used ini- 
tial fluctuations in the density and temperature of 
rms amplitude equal to 0.1 times their mean values 
(=1), with zero initial velocity fluctuations. For the 
full ISM-like runs, we have used initial fluctuations 
with rms amplitudes equal to the mean values of the 
variables, except for the magnetic fleld, whose fluctu- 
ations are taken as described in 

3. Results and discussion 

3.1. Development of the pure instability 

3.1.1. Morphology and dynamics 

The density fields for runs PFS and PFS at an ad- 
vanced stage of the instability development are shown 
in Figs. |l|a (left panel) and |l]b (right panel), respec- 
tively. In run PF5 (0a), the gas is seen to evolve into 
a "bee-nest" structure formed by a network of dense 
(~ 10 cm"'^), cold (~ 50 K) filaments enclosing warm 
(- IQ-* K) and rarefied (~ 0.3 cm'^) cells. The cells' 
shape and size are determined by the initial condi- 
tions. In this run we have taken = 2, implying 
that the characteristic scale of the initial fiuctuations 
is one half of the box size. Note that the TI leads 
to the collapse of regions containing a density excess, 
whose characteristic size is one fourth of the integra- 
tion domain. The roundish blobs formed at intersec- 
tions between filaments originate from density max- 
ima, while the dense filaments can be traced back to 
"saddle points" located between two density maxima 
and two minima in the initial conditions The evo- 
lution of this simulation is reminiscent of "pancake" 
formation due to the gravitational instability in cos- 
mological simulations of large-scale structure forma- 
tion (e.g., Padmanabhan 1993). However, this bee- 

^The initial conditions for the 2D density field are mosaics of 
small-amplitude, alternating maxima and minima. Typically, 
imaginary lines joining two maxima and lines joining two min- 
ima intersect, forming an "x" pattern, at the center of which 
there is a saddle point. 



nest structure is transient, as we discuss below. 

The condensation process in principle stops after 
the pressures in the warm and the dense (cold) phases 
equalize and any remaining kinetic energy due to the 
inertia of the gas is dissipated. This implies that the 
mass redistribution induced by the instability does 
not end when the density has "crossed" the transition 
values, but at more disparate values of the density 
such that the pressures in the two phases are equal. 
For example, fig.j3 shows the equilibrium pressure, 
as given by eq. (Q) for runs PF5 and Pr3. It can 
be seen that in order for the pressure in the warm 
phase to equal that of the dense phase, a density of 
roughly 0.03 cm~^ is required in the former. (Note 
that, since run PF5 has 712 = 0, the equilibrium can 
only be reached by decreasing the warm phase' den- 
sity.) However, at the final time reached by run PFS, 
the pressure in the warm phase is still much larger 
than in the dense phase, as can be seen in fig. ^. That 
is, this run stopped only because the shocks became 
too strong for the code to handle, but it was still in a 
highly dynamical contracting stage. 

The highly dynamical nature of the condensation 
process is clearly seen in fig. ^ which shows horizon- 
tal cuts at ?/ = 42 through the (log of the) density 
(dotted line), (log of the) pressure (dashed line), and 
the a;- velocity (solid line) fields of run PF5 (fig. |l|a). 
The cuts go through two blobs and one filament. A 
double-shock structure is observed in the blobs, which 
arc the two-dimensional equivalent of one-dimensional 
shock-bounded slabs. This should also be the case of 
the filament, but it is too thin for the code to resolve 
the double-shock structure, since shocks in the code 
are spread over roughly five grid points by viscosity. 

In order to see the subsequent evolution, we thus 
performed another run, named PF.25 with slower 
heating and cooling rates (by decreasing the coeffi- 
cients) but with the same temperature dependence 
as run PF3 (except for a slight change in T2), which 
therefore evolved more slowly, developing milder shocks, 
and allowing us to reach more advanced stages of the 
instability development. A temporal sequence in the 
evolution of this run is shown in fig. |^. It is seen that 
this run also develops the bee-nest structure, but later 
the filaments disrupt and are accreted into the actual 
peaks. Moreover, some of the cells "implode", and 
their surrounding filaments and peaks collapse into 
the void, due to the small pressure differences be- 
tween cells, traceable to the randomness of the initial 
conditions. The behavior thus continues to resemble 
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that induced by gravitational instability in cosmolog- 
ical simulations, as the pancakes later collapse into 
isolated structures (e.g., Padmanabhan 1990). 

Run Pr3 (fig. |l|b), on the other hand, does not 
show the bee- nest structure of run FFS, but shows 
roundish blobs within elongated, isolated structures, 
and only one truly filamentary structure. More im- 
portantly, run PrS never develops the thin, dense fila- 
ments seen in run PT5. Instead, only thick, moderate- 
density bands appear between the peaks, and then 
disappear. We attribute this difference in part to the 
fact that in run PFS, p3 is very close to the mean 
density (see Table 1). This causes any initial under- 
density to quickly reach the stable range, implying 
that a relatively large amount of mass remains in the 
warm phase (note that the buildup of an underdensity 
entails the evacuation of that region) . 

However, the main reason for the non-appearance 
of the filaments seems to be again the fact that the 
densities of the warm and cold phases required for 
equilibrium are significantly more distant from the 
mean density than p2 and p^. This implies that, dur- 
ing the development of the instability, the pressure 
maxima do not occur at the centers of the voids, but 
in the boundaries between the voids and the high den- 
sity regions (fig. ^). In the particular case of run PFS, 
P3 is so close to (p) that the voids reach smaller pres- 
sures earlier than the unstable filaments. The latter 
then reverse their condensation process, loosing their 
mass to both the voids and the peaks, and ultimately 
merging with the intcrcloud medium. The peaks do 
not suffer this "rebound" because, due to their larger 
initial densities, they reach lower pressures faster than 
the voids. 

Even though the above effect is particularly no- 
ticeable in run PPS, we expect it to apply in general. 
Therefore, unless the density field is perfectly uniform 
within the phases, and the interface between them is 
a perfect discontinuity, a larger pressure is expected 
to be always present near the interface between the 
two phases than in either one of them. This suggests 
that a steady state with constant pressure everywhere 
is difficult to attain. 

3.1.2. On the final state of the instability 

Another consequence of the above discussion is 
that, since the gas is a fluid, the boundary of the 
"clouds" that make up the dens e phase are accretion 
shocks, as discussed by BVS99, rather than the qui- 



escent density jump at constant pressure usually as- 
sumed in the hterature. Of course, quiescent pressure 
equilibrium may be established after the pressures 
have equalized, and the boundary shocks subside, but 
this may require very long times, and besides, the bal- 
ance is very delicate. In general, there is a continuum 
of possible equilibria between the dense and rarified 
phases, obtained by equating their pressures, as given 
by eq. (||), and solving for pd, the density in the dense 
phase, as a function of p-c, the density in the rarified 
(warm) phase. We obtain 



Cl2 



l//3i 



l/fe 



^ 1/71 



(12) 



The dependence of p^ on p-^ is thus determined essen- 
tially by the ratio 734/712. As an example, plugging 
in the corresponding parameters for run PFS, which 
has 712 = 1/3, we find that pd = 3.2 x Iffip^'^^. 
Thus, in this case, the equilibrium density of the 
dense phase depends sensitively on the density of the 
warm phase, minute fluctuations in pr requiring large 
changes in p^. For example, at pr = 0.1, pd = 16.8, 
while at p-^ — 0.3, pd = 205 (all in units of the mean 
density). Such large density variations for the dense 
phase's density will require highly dynamical adjust- 
ments which will most likely involve new shocks, espe- 
cially since the sound speed in the warm phase is much 
larger than in the dense medium. This, together with 
the result by Murray et al. (1993) that the condensed 
clouds are likely to be set in motion relative to the 
diffuse phase by buoyancy or other effects, and then 
easily disrupted by Kelvin-Helmholtz-type instabili- 
ties, strongly suggests that the static conflguration is 
highly improbable. 

In cases where 712 = 0, eq. ( [T^ ) does not apply, and 
the pressure of the dense phase becomes independent 
of its density. Thus, pressure equalization can only be 
accomplished through evacuation of the warm phase 
until its pressure drops to the value within the dense 
gas. In the meantime, the density of the dense gas 
may increase by very large amounts, most likely in- 
volving violent compressions. We conclude that the 
static configuration is very hard to realize, in agree- 



ment with the conclusions of BVS99 



In any case, the formation of static configurations 
is possible in principle, although we have seen here 
that this is a direct consequence of the change in the 
effective polytropic exponent 7° from one phase to 
the other, as discussed in detail by BVS99| . This sug- 
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gests that cores within molecular clouds, which have 
essentially the same value of 7** as their parent clouds, 
should not be expected to reach hydrostatic equilib- 



3.1.3. PDFS 

The density PDF has been the subject of much re- 
cent work. Vazquez-Semadeni (1994) and Padoan, 
Nordlund & Jones (1997) have reported lognormal 
PDFs for two- and three-dimensional isothermal flows, 
respectively. Passot & Vazquez-Semadeni (1998) dis- 
covered that the functional form of the PDF for poly- 
tropic {P oc p'^ ) one-dimensional flows depends on 
the effective polytropic exponent 7°, developing a 
power-law tail at high densities for < 7*^ < 1 and 
viceversa, with a lognormal indeed developing when 
7'' = 1 (see also Nordlund & Padoan 1999). Scalo 
et al. (1998) found that the simulations of the ISM 
by Passot et al. (1995), as well as 2D Burgers-like 
(without thermal pressure) runs developed power-law 
PDFs analogous to the < 7° < 1 cases of Passot & 
Vazquez-Semadeni (1998). 

In the thermally unstable case, the instability causes 
the evacuation of material from the unstable to the 
stable regions, producing a two-phase density field 
characterized in principle by a bimodal PDF. In the 
remainder of this paper we use the density PDF as a 
diagnostic of the relative importance of the instability 
in the dynamics of the various simulations. 

The density PDFs for runs PT5 and PFS corre- 
sponding to the fields shown in fig. |l| are presented in 
fig. ^ The vertical lines denote the transition densi- 
ties p2 and p3 for the two runs. The morphological 
and dynamical features described in the last section 
manifest themselves in the PDFs as well. In both 
cases, a peak in the PDF is observed just below the 
lower transition density, p^. However, in the case of 
run Pr5, p3 is not too similar to the mean density, and 
this implies that a substantial amount of mass has to 
be transferred to the dense phase. As a consequence, 
a noticeable peak appears in the PDF largewards of 
the upper transition density, p2. Instead, for run PFS, 
much of the mass remains in the low density phase, 
and there is no clear peak above p2. Only a change 
in the slope of the PDF is seen there. 

Furthermore, the PDF for run PFS seems to extend 
to much higher densities than that of run PFS. We 
interpret this as a consequence that run PFS reached 
a more advanced stage in the development of the in- 



stability than run PFS, so that densities much higher 
than the transition value are reached as the densities 
in each phase approach the equilibrium values. In 
fact, this phenomenon should have also been observed 
in run PFS, had it reached more evolved stages. In- 
deed, in fig. [s] we show the density PDF at three differ- 
ent times for run PF.25 which, as mentioned above, is 
similar to run PFS but evolves more slowly, producing 
weaker shocks and therefore reaching more advanced 
stages before the code cannot handle the shocks any- 
more. In this figure it is clearly seen that the PDF 
is similar to that of run PFS at first, but then over- 
shoots and becomes more similar to that of run PFS, 
the peak above the upper transition density disap- 
pearing altogether. 

We conclude from this section that both a bimodal 
PDF and a unimodal one with a slope change above 
P2 may be considered the signatures of the TI. 

3.2. Effects of turbulent forcing 

We now discuss the effects of adding turbulent forc- 
ing to purely hydrodynamic simulations in the pres- 
ence of the TI. Our main objective is to determine 
whether a turbulent regime can prevent the develop- 
ment of the instability, in the sense of destroying the 
bimodality (or slope change; cf. §3.1.3) of the density 
PDF, so that the two-phase nature of the medium 
disappears. In particular, we also wish to determine 
whether different types of forcing have different effects 
on the development of the instability. 

To this end we consider run PFS, which is only 
subject to the TI, and three forced runs with ex- 
actly the same heating and cooling functions, but sub- 
ject also to various types of forcing. As described in 
sec. |2.3| , we refer to the forced runs generically by 
"H" , and then add a mnemonic for the type of forc- 
ing used. Run HPS. 5 uses large-scale random forc- 
ing with S0% solenoidal (rotational) content ("PS" 
— "percent solenoidal"). Run HPSl uses the same 
type of forcing, except with 100% solenoidal content. 
Finally, run HSF ("star formation") uses small-scale 
(a few pixels) stellar-like forcing, which is 100% com- 
pressible. 

In fig. II we compare the density PDFs for the four 
runs. Interestingly, the two runs with large-scale forc- 
ing are not able to counteract the development of the 
instability, as indicated by the fact that their PDFs 
still show a significant slope change above the up- 
per transition density p2, similar to that of run PFS 
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(§3 1.3|). In fact, in these cases the turbulence appears 



to reinforce the instabihty, since both runs stopped 
due to the large gradients at earlier times that the 
pure-instability run PFS. Note, however, that a peak 
above p2 is not present for either run, apparently be- 
cause the turbulence contributes to the overshooting 
effect, which spreads the mass in the dense phase over 
a larger range of densities than in run PFS. The in- 
ability of the turbulence to counteract the instability 
manifests itself also in the density fields of both HPS 
runs (not shown), which show essentially the same 
type of filament network as run PFS, but distorted 
by the large-scale turbulent motions. 

Run HSF, on the other hand, due to the small-scale 
and "explosive" nature of the forcing used, is able to 
counteract the collapse of the filaments, since the stel- 
lar heating reverts the compression of the dense gas 
by the warm phase. This effect causes the PDF for 
this run (dash-dot) to exhibit the interesting feature 
that the high-density peak spreads out, and in partic- 
ular invades the unstable density region line in fig. ^, 
ultimately erasing the signature of the TI. Indeed, in 
this case, the shoving of the gas by expanding shells 
continually redistributes it among the various thermal 
regimes. This erasure effect is clearly shown in fig. |l^, 
in which the run's PDF is seen to transit from clear 
bimodality to essentially a power-law tail at densi- 
ties above the mean. However, this run is "unstable" 
in the sense that it displays a runaway star forma- 
tion rate (recall star formation in the simulations is 
triggered by large densities). This suggests that TI 
in this run still has the effect of promoting excessive 
star formation. 

The inability of large-scale forcing to prevent the 
development of the instability may be understood in 
terms of the time scales associated with the insta- 
bility's growth rate and with the nonlinear turbulent 
crossing time. In the absence of self-gravity {J — 0), 
magnetic field {Bo = 0) and rotation (il = 0), the 
growth rates given by eqs. (|l0|) and ( ITTI ) become iden- 
tical, giving uj"^ — — (7°i i+i/7)7'cqfc^ = 0.41fc^, using 
Tcq = 0.23, 7°i^i+i — —3 and 7 = 5/3; i.e., the growth 
rate increases linearly with wavenumber. This im- 
plies that the characteristic time scale for the growth 
of fluctuations under the effect of the instability tti 
increases linearly with the size scale I. 

On the other hand, for a turbulent flow, the veloc- 
ity difference across a scale I varies as oc where 
a = 1/3 for a Kolmogorov spectrum, and a = 1/2 for 
a shock-dominated flow (see, e.g., Vazquez-Semadeni 



et al. 2000 ). This implies that the nonlinear crossing 
time at scale I varies as tnl oc Z'^, with /3 — 2/3 
for a Kolmogorov flow, and /3 — 1/2 for a shock- 
dominated flow. In either case, it is seen that, for 
sufficiently small scales, tti < tnl- Thus, a turbulent 
flow forced at large scales, and cascading its energy to 
small scales, cannot prevent the development of the 
instability at sufficiently small scales. 

Note that the above argument assumes that the 
effect of turbulence is to always fight the instability. 
However, for compressible flows, it is in fact possi- 
ble that in a fraction of the volume the turbulent 
transfer from large to small scales occurs via the for- 
mation of shock-bounded slabs (e.g., Elmegreen 1993; 
Vazquez-Semadeni et al. 1995| ; |BVS99| ), in which case 



the transfer from large scales promotes compression at 
small scales, thus helping the instability rather than 
opposing it. In fact, Hennebelle & Perault (1999) 
have recently shown the triggering of TI by compres- 
sive motions in originally stable flows. In either case, 
it is seen that in a turbulent flow forced at large scales, 
the turbulent transfer is incapable of preventing the 
development of the instability in general. It is worth 
noting that these results are qualitatively similar to 
the inability of large-scale turbulent forcing to sup- 
press the gravitational instability (Leorat, Passot & 
Pouquet 1990). On the other hand, as we have seen, 
stellar-like small-scale forcing, applied at the sites of 
maximum density, is capable of erasing the signature 
of the instability. 

3.3. ISM-like simulations 

We now consider a set of three ISM-like turbulent 
simulations including self-gravity, the magnetic field, 
large-scale shear, the Coriolis force and energy input 
due to star formation, but differing in the cooling 
functions in order to make them either thermally sta- 
ble, marginal or unstable. These runs are respectively 
labeled ISMS, ISMM and ISMU. The initial condi- 
tions include also velocity and magnetic field fluctua- 
tions. In particular, run ISMU has identical heating 
and cooling functions as run PF3. 

It should be pointed out that we have tested nu- 
merically that run ISMU is also unstable with re- 
spect to the combined effect of all physical agents 
inc luded, but withou t sta r formation. As di scussed 
by Elmegreen (1994) and Passot et al. 1995 , in the 
presence of shear, eqs. (|lO|) and (|ll|) are only valid 
near t = 0. At later times, the development of the 
perturbations has to be investigated numerically. We 
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also tested that, in the absence of shear, run ISMU is 
unstable as well, as indicated by the dispersion rela- 
tions. 

Figure |l^ shows the density PDFs for these three 
runs. These runs reach a stationary regime, and 
therefore we have integrated the PDFs over several 
timesteps to improve the statistics. The most impor- 
tant difference in the PDFs actually occurs at low 
densities. This seems to be a consequence that the 
lower transition density for run ISMU is very close 
to the mean density, as was the case for run PFS, 
and thus the warm phase is not strongly evacuated. 
Instead, for the marginal and stable runs, the lower 
transition densities are quite small (see Table 1). Al- 
though at first one might think the transition den- 
sities should not matter in non-unstable cases, they 
actually do. This is because the density range be- 
tween p2 and P3 continues to be "softer" (i.e., more 
compressible) than densities below ^3. Thus, once the 
density drops to these values, the thermal pressure 
acts more strongly against further evacuation of these 
regions by the turbulence (cf. Passot & Vazquez- 
Semadeni 1998). 

At large densities, the three curves are remarkably 
similar, with only a trace of a slope change in the PDF 
of the unstable run above its upper transition density. 
However, this slope variation is minimal compared to 
the dramatic one seen in the PDF of run PFS (fig. 

This shows that the effect of TI on the density 
distribution has been erased by the small-scale energy 
injection from the stellar sources. 

4. Summary and conclusions 

In this paper we have discussed the development 
of the thermal instability (TI), both by itself and in 
the presence of turbulent forcing. We have also ex- 
amined its role in determining the statistical distri- 
bution of the density field, as measured by its proba- 
bility density function (PDF), in an ISM-like regime 
with magnetic fields, self-gravity, stellar-like energy 
injection, and rotation. We used two-dimensional nu- 
merical simulations on the plane of the Galactic disk 
at moderate resolutions, allowing us to perform a rea- 
sonable coverage of parameter space. Our simulations 
are not able to reach the stage at which pressure equi- 
librium is established because of the strong shocks 
involved, but stop instead at some earlier time. 

The simulations of the development of the TI 
alone, starting from small perturbations, showed that 



the morphologies that arise depend quite sensitively 
on the proximity of the transition densities p2 and 
Pa (the values of the density bounding the unstable 
regime from above and below, respectively) to the 
mean density (p) of the medium. If the lower tran- 
sition density is very close to the mean density, a 
large fraction of the mass remains in the warm (low- 
density) phase, and there is little mass available for 
forming the dense phase ("clouds"), which consists es- 
sentially of isolated, roundish structures. Conversely, 
simulations in which pa is significantly smaller than 
(p) evacuate more mass from the warm phase, and 
develop a transient filament network before the fila- 
ments are accreted into the actual peaks, similarly to 
pancake formation in cosmological large-scale struc- 
ture formation simulations. The PDFs of advanced 
stages of the TI are either bimodal with peaks above 
and below the upper and lower transition densities, 
respectively, or else show the low-density peak plus a 
slope change above p2. 

The condensation process which originates the "clouds' 
was shown to be highly dynamical, with their bound- 
aries being accretion shocks rather than static density 
discontinuities. Although in principle the formation 
of the latter is possible, it was shown that, for real- 
istic cooling functions, the equilibrium density in the 
cold phase depends sensitively on that of the warm 
gas, so that small changes in the latter require large 
changes in the former, which most likely would imply 
the formation of shocks. Furthermore, if the density 
is continuous from one phase to the other, then the 
two phases must be mediated by larger-pressure re- 
gions, which are expected to induce further motions. 
The static configuration is possible in principle, but 
appears highly unlikely. 

The inclusion of energy- injection processes ("forc- 
ing") has very different effects depending on the na- 
ture of the forcing. Random, large-scale forcing is not 
able to inhibit the development of the instability, and 
only distorts the structures it forms on time scales 
longer than the growth rates of the instability. In 
fact, large-scale turbulent modes naturally generate 
compressions, so the large-scale forcing actually aids 
the instability. Instead, small-scale, stellar-like forc- 
ing systematically injects energy within the structures 
formed by the instability, and is capable of destroy- 
ing the signature of the instability in the PDF, push- 
ing gas from the dense phase back into the unstable 
regime. However, the flow continues to be extremely 
compressible, as indicated by an ever-growing star- 
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formation rate and the inability to develop a station- 
ary regime. 

Finally, ISM-like simulations with and without TI 
show little difference in their PDFs, suggesting that 
the combined effect of the stellar-like forcing, the 
magnetic pressure and the Coriolis force overwhelm 
the thermal pressure deficit in the unstable cases. 
Furthermore, all of the ISM cases reach stationary 
regimes, indicating that the presence of a TI in the 
medium is of relatively minor importance in the over- 
all cycle of such a regime. However, since the turbu- 
lent forcing originates from the stellar activity, the 
importance of the TI in the formation of the first 
generation of stars in a galaxy cannot be ruled out, 
although, on the other hand, the weaker cooling in 
those epochs due to the absence of heavier elements 
may have prevented the appearance of the instability 
altogether. The likely disruption of clouds by Kelvin- 
Helmholtz instability (Murray et al. 1993) suggests 
that TI, even if it can occur in primordial galaxies, 
will be unable to form stable clouds. 

We conclude by emphasizing that a bimodal tem- 
perature is not necessarily evidence of a multi-phase 
medium, but only of a nonlinear dependence of the 
equilibrium temperature on density. A true multi- 
phase medium must involve clearly distinct phases, 
which, in the case of TI, means different thermal equi- 
libria, some of which may be stable and others unsta- 
ble. For example, Korpi et al. (1999) have also inves- 
tigated the temperature and density PDFs in realistic 
(but not self-gravitating) 3D ISM simulations, but in 
their model there is no background heating, so in ef- 
fect their simulations are always in a cooling regime 
(albeit possibly very slowly at low densities), except 
in the vicinity of stellar thermal energy sources, and 
thus there are no stable "phases" . Yet, those authors 
still interpret the resulting bimodal temperature dis- 
tribution as an indication that the medium is "mul- 
tiphase" , in spite of the fact that there are no stable 
phases in their model, the pressure varies by over one 
order of magnitude, and the density pdf shows no sign 
of bimodality. A bimodal temperature distribution is a 
simple consequence of the form of the cooling law, and 
will occur in any model of the ISM in which gas can 
attain large temperatures and cool; even in a nondy- 
namical time-dependent model with no isobaric insta- 
bilities (see Gerola, Kafatos, & McCray 1974). But, 
in general, this does not imply the existence of sta- 
ble or unstable phases in the medium. Thus, we feel 
that the concept of a "multiphase" medium has some- 



times been inappropriately used in the ISM literature. 
For example, we see no rationale for referring to the 
correlation function turbulence model of Norman & 
Ferrara (1996) as a "continuum of phases", when in 
actuality it is simply a non-isothermal density con- 
tinuum, while, on the other hand, such terminology 
perpetuates the idea that the ISM is really controlled 
by phase transitions. 
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"Turbulencia Interestelar" to EV-S and NASA grant 
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Fig. 1. — Density field of runs a) PFS {left) and b) 
Pr3 (right). In tliese runs, only the TI is at work, 
starting from small density fluctuations. Run PFS is 
seen to develop a network of filaments, while run PFS 
transits directly to the formation of isolated density 
structures. 




Fig. 3. — A pixel-by-pixel plot of pressure vs. density 
for run PF5. Compare to the solid curve in fig. |^. The 
thickness of the negative-slope portion is due to the 
finite cooling times, which cause the temperature of 
a fiuid parcel to depend not only on its instantaneous 
density, but also on the speed at which it is varying. 
Numerical noise due to the abrupt transition from 
one cooling regime to another also contributes to the 
scatter. 
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Fig. 2. — Equilibrium pressure as a function of den- 
sity as given by eq. (^, for the cooling curves used in 
runs PF5 [solid line) and PF3 [dotted line). 
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Fig. 4. — Horizontal cuts through run PFS at y = 42 
of the density (dotted line), x- velocity (solid line) and 
thermal pressure (dashed line). The two structures 
on the left are blobs (see fig. and the one on the 
right is a filament. Double-shock structure character- 
istic of shock-bounded slabs, with a central plateau 
bounded by abrupt jumps, is seen in all three vari- 
ables in the blobs. It cannot be seen in the filament 
on the right because this structure is too small, and 
viscosity blends the two shocks in a single one. 
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Fig. 5. — Three snapshots of the density field of run 
Pr.25, at t = 12.4 (top), t = 16.4 (middle) and 
t = 19.3 (bottom). The times are in code units, with 
1 code time unit — 13 Myr. Note that the filament 
network fragments and starts producing isolated en- 
tities. 
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Fig. 6. — Images of the density {left) and thermal 
pressure (right) for run PT3 at an intermediate stage 
of the instabihty development. White means large 
values and black means low values. Although the 
voids have larger pressures than the density peaks, the 
largest pressures occur at the interfaces between voids 
and filaments or peaks. Filaments have larger pres- 
sures than the voids, and therefore "dissolve" rapidly. 
The presence of such high-pressure interfaces is ex- 
pected in general in the development of the instability 
if the density field is not perfectly smooth. 
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Fig. 7. — Probability density functions (PDFs) or his- 
tograms of the density field for runs PFS {solid line) 
and PF3 {dotted line). The transition densities for 
both runs are shown as vertical lines, with the same 
line type as the histogram for each run. The PDF 
of run Pr5 shows two peaks outside its correspond- 
ing transition densities, while that of run PFS shows 
only the low-density peak. Instead of the high-density 
peak, a slope change beyond its upper transition den- 
sity is seen. We consider either kind of PDF a signa- 
tured of the dynamical effect of the TI. 
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Fig. 8.— Evolution of the density PDF of run Pr.25. 
The PDFs are shown at the same times as the fields 
shown in fig. ^. The PDF is initially bimodal, but, as 
the instability continues to develop, the high-density 
peak is stretched out towards even higher densities, 
and becomes simply a flatter portion of the PDF. 



Fig. 9. — PDFs for four unstable runs with exactly 
the same cooling functions, but with different energy- 
injection (forcing) regimes. The transition densities 
are indicated by the vertical lines. The unforced run 
is run PPS, adopted as a reference. Run HPSl has 
large-scale, purely solenoidal forcing; run HPS. 5 has 
large-scale forcing with 50% solenoidal content; fi- 
nally, run HSF has small-scale, stellar-like forcing. 
Only the latter run is able to destroy the signature 
of the instability. 
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Fig. 10.— Evolution of the PDF of run HSF. Be- 
fore star formation (SF) starts, the run exhibits the 
normal bimodal PDF. As SF injects energy, the high 
density peak is spread both to high and low densities, 
until finally a single power-law PDF remains at high 
densities, destroying the signature of the instability. 
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Fig. 11. — PDFs for three ISM- like simulations in- 
cluding stellar energy injection, the magnetic field, 
the Coriolis force and self-gravity. The PDFs are av- 
eraged over several timesteps to improve the statis- 
tics. The runs differ in the cooling functions, which 
render them either stable, marginal or unstable. The 
transition densities for the unstable run are indicated 
by the vertical lines. The main difference between the 
PDFs is at low densities, due to the fact that its lower 
transition density is higher than those of the the other 
two runs. At high densities, virtually no difference is 
seen between the PDFs. 
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Table 1 
Run Parameters 



Run # 


Label 












£23,34^ 


Tag 




P3' 


T j 

eg 




7I3' 


Forcing™ 


ISM" 


1 


pro5 


5 


1 


.25 


4.18 


80.62 


7.184 


.0398 


1.56 


.696 


.2347 


0.0 


-3 






2 




3 


1.5 


.25 


4.18 


170 


3.022 


.0398 


2.22 


.993 


.9711 


.33 


-3 






3 


pro. 25 


.25 


1.5 


.25 


4.18 


200 


.6325 


.01 


1.25 


.395 


.0244 


.33 


-3 






4 


HPS. 5 


5 


1 


.25 


4.18 


80.62 


7.184 


.0398 


1.56 


.696 


.2347 


0.0 


-3 


ps = .5 




5 


HPSl 


5 


1 


.25 


4.18 


80.62 


7.184 


.0398 


1.56 


.696 


.2347 


0.0 


-3 


ps = I 




6 


HSF 


5 


1 


.25 


4.18 


80.62 


7.184 


.0398 


1.56 


.696 


.2347 


0.0 


-3 


SF 




7 


ISMS 


3 


1.5 


1.5 


3.0 


150 


150 


.0398 


2.52 


.020 


.0737 


.33 


.33 


SF 


on 


8 


ISMM 


3 


1.5 


1.0 


3.0 


160 


31.91 


.0398 


2.36 


.094 


.0940 


.33 


0.0 


SF 


on 


9 


ISMU 


3 


1.5 


.25 


4.18 


170 


3.022 


.0398 


2.22 


.993 


.9711 


.33 


-3 


SF 


on 



^Diffuse heating rate 

''Cooling law exponent between 10 K and T2 
'^Cooling law exponent between T2 and 10'* K 
'^Cooling law exponent between 10* K and 10^ K 
"Cooling law coefficient between 10 K and T2 

^Cooling law coefficients between T2 and 10* K and between 10* K and 10^ K 

^Transition temperature in units of 10* K 

'^Transition density at T = T2 

'Transition density at T = 10* K 

j Equilibrium temperature at (p) in units of 10* K 

'^Effective polytropic index between 10 K and T2 

'Effective polytropic index between T2 and 10* K 

™Type of forcing 

"ISM additional physics (magnetic field, self-gravity and rotation) 
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